The role of the diffusion in the predictions of the classical nucleation theory for quasi-real systems differ in dipole moment value

In this paper, we examine the crystallization tendency for two quasi-real systems, which differ exclusively in the dipole moment's value. The main advantage of the studied system is the fact that despite that their structures are entirely identical, they exhibit different physical properties. Hence, the results obtained for one of the proposed model systems cannot be scaled to reproduce the results for another corresponding system, as it can be done for simple model systems, where structural differences are modeled by the different parameters of the intermolecular interactions. Our results show that both examined systems exhibit similar stability behavior below the melting temperature. This finding is contrary to the predictions of the classical nucleation theory, which suggests a significantly higher crystallization tendency for a more polar system. Our studies indicate that the noted discrepancies are caused by the kinetic aspect of the classical nucleation theory, which overestimates the role of diffusion in the nucleation process.


Results
Similarly, to our previous examinations of the system I, we began the studies of system II from the constructing the perfect FCC crystal structure constructed from 2048 molecules and heating it from 10 K up to the temperature which is about 50 K higher than the temperature at which we observe a rapid increase in the volume, see Fig. 1.
On the basis of our recent results for the system I, we can state that the rapid increase in the volume from about 0.086 to 0.094nm 3 , which is observed around T = 220K , indicates on the melting of the crystal structure. Therefore, the thermodynamic conditions at which the system II is in the liquid phase can be recognized in the www.nature.com/scientificreports/ similar way, see increase in volume from about 0.088 to 0.101nm 3 at T = 420K. Next, we cooled the systems up to the starting temperatures. During this process the drop in the volume can be observed. The latter indicates on occurrence of the crystallization process, which takes place at T cr = 130 and 280K , respectively for the systems I and II. According to CNT the nucleation rate is expressed as follows 13 where ρ liq is the number density of the liquid, D is a diffusion constant, k B is the Boltzmann constant, and �W = 16 3 γ 3 (�Gυ ) 2 is the nucleation barrier, in which �G υ denotes the driving force per volume unit (i.e., the difference between Gibbs free energy for liquid and bulk phases) and γ is the interfacial free energy (IFE). The next physical quantity determining the occurrence of the overall crystallization process is the crystal growth rate, role of which can be computed by the following expression where A U (T) describes the molecular mobility and can be approximated by D · a/ 2 , in which a is the average width of the crystal lattice spacing ( a ≈ ρ 1 3 cr , ρ cr denotes the number density of the crystal), ≈ ρ −1/3 liq is the atomic jump distance, whilst f (T) describes the grow mechanism, which for the normal growth ≈ 1. For all thermodynamic conditions at which during the cooling systems remain in the liquid phase the diffusion constant is determined using the GROMACS software on the basis of the mean square displacement calculated for atoms for long times. In this way both translational and rotational contributions to the molecular motion are considered. Subsequently, the estimated dependences D(T) are approximated by the Vogel-Fulcher-Tammann equation for the needs of further analysis (see Fig. 5a, which we discussed later). To estimate the value of G we use the method proposed by Gutzow 52 with the redefined integration pathways 53 , according to which, at isobaric conditions, the driving force for the crystallization takes the following form �G(T) = − T T m �S(T)dT , where S is the difference in the entropy between the liquid and the solid phases. Taking into account the melting boundary conditions �S m = �H m /T m , S can be calculated using obtained directly from simulation-runs values of the enthalpy H(T) and the classical relation between enthalpy and the entropy, T = (∂H/∂S) p . The estimation of the melting temperature has been done using the liquid-solid coexistence method. In this order, we visualized the structure up to which each system crystallizes. Then, we determined the fragments characterized by a high degree of order, which for both systems are characterized by the triclinic shape and consist of molecules placed in corners. Based on the latter, we constructed another crystal structure and equilibrate it at the temperature close to T cr for both systems. Since we observed that the small defects occur again, we selected the set of 5 × 5 × 5 molecules within which the created crystal structures were highly ordered and those crystal fragments are used for further examination. On their basis, we construct the crystal structures consisted of 2250 molecules, and equilibrate it at the temperatures significantly smaller than T cr . Subsequently, we heat the systems to confirm that the crystals are stable at higher temperatures. The results are presented in Fig. 1 of the main text. It can be seen that created crystal structures do not tend to melt, although, during the heating process, the tiny and step increase in volume is detected for both systems. Probably, these changes of volume are results of the transformation to the different polymorphic form, which is more stable at higher temperatures. This scenario seems to be supported by the evident visible change in the temperature dependence of volume, which is observed around www.nature.com/scientificreports/ 100K during the cooling of the system II. However, to confirm this suspicion further researches are required. Nevertheless, it is worth mentioning that created structures are more stable at higher temperatures than the ones resulted from starting FCC configuration. This fact encourages that the crystal structures established by us can be used to determine T m . Hence, we constructed the special biphasic simulation box, within which 3456 molecules had been equally divided between the crystal and liquid phases separated by a small gap. Due to the fact that at melting conditions, the solid and liquid phases remain in the thermodynamic equilibrium, T m can be determined by the examination of the behavior of the biphasic system. The performed simulations of biphasic box last for 5ns and have been done at range of temperatures differ by 1K . Following the visual examination of the obtained configurations we determine that melting temperatures, which equal T m = 150K and T m = 326K for the systems I and II, respectively. At this point, we have to comment that in our previous studies 51 we determined the T m for the system I, using liquid-solid coexistence method as well, and we obtained that T m = 194K . However, in that experiment we did not determine the crystal structure. Instead of that, we employed the structure up to which the liquid had crystalized. Consequently, we probably employed the structure stable at higher temperatures instead of the one, which is the most energetically optimal. However, in this work, we intend to focus on the most fundamental case, i.e., we estimate the crystallization tendency against the desired (and the most energetically optimal) structure. In this context, it is worth noting that the values of T m estimated herein are in similar relations to those of T cr and also to the temperature at which initial crystals melt, which suggests that the determined structures are mutually appropriate. However, the most challenging is the determination of the γ . Fortunately, the special computational method for calculation of γ have been proposed. The two main approaches are the cleaving potential method 54-56 and the capillary fluctuation method [57][58][59][60][61] . Despite that both methods are applicable only at the melting conditions, they strongly differ in the way of work. In the cleaving potential method, the biphasic solid-liquid system is transformed into two separate systems (liquid and solid) by means of external potentials. Then, γ is estimated on the basis of the work which is performed by those potentials during the transformation process. However, the precise application of this method is associated with some technical difficulties. It is because the reversibility of the transformation process must be ensured, and therefore, the accurate control on the transformation process is needed 62 . Alternatively, γ can be calculated in the more direct way using the capillary fluctuation method (CFM), which requires only one simulation run, during which any knowledge of the complex process of the interface creation from separated bulk systems is not needed. Instead of that, through the simulation run, the fluctuations of the interface are measured. It enables the estimation of the interface stiffness, which is related to γ . The remarkable advantage of CFM is the fact that it considers the anisotropy of the interface, whereas the cleaving method is recognized as more accurate. Till now, both methods have been applied to calculate γ values for model systems such as hard-spheres 55,63 and Lennard-Jones 25,56,60 . It must be however noted that for the real materials the CFM is more often employed, which is mainly stimulated by the ease of its application. Consequently, using the CFM the γ values have been calculated for metallic compounds [57][58][59]64 , alloys 65,66 , and a few molecular systems 67 including pharmaceuticals [68][69][70] . Hence we decided to employ CFM to determine γ for studied herein systems. The use of CFM requires creation of the biphasic box. However, the considered solid-liquid interface must be the quasi-one-dimensional, and therefore the special geometrical conditions of the simulation box have to be ensured, i.e., when interface is perpendicular to the length of the system, L x , its thickness must be much smaller than its width, L z ≪ L y . Then the interface fluctuates only in the one dimension ( x ). Consequently, we construct the box containing 5000 RLM divided equally to the crystal and liquid phases. It is also worth mentioning that due to boundary conditions the simulation of the biphasic system implies the existence of two interfaces, of which the fluctuations magnitudes are studied. The convenient way used to determine the temporary position of the interface is the calculation of the rotational-invariant order parameter 64,71-75 ( RIOP ) for geometrical center of the molecules. The RIOP enables the distinction between solid-like and liquid-like molecules, because the solid-like molecules are characterized by the significantly higher values of the order parameter. The example of obtained results is presented in Fig. 2a, where the calculated RIOP for each molecule of system I is plotted as a function of the position of the molecule in the dimension perpendicular to the interface plane ( x direction).
As we already mentioned the liquid-like and the solid like molecules can be clearly distinct. The evolution of the RIOP can be described by the following function RIOP y = o s +o where o s,l are the average values of the RIOP in the solid and liquid, δ 1,2 are effective widths of the interfaces, and h 1,2 (x) are functions describing the positions of the interfaces in capillaries, i.e., sections from x to x + x , which are orthogonal to the interface. During the simulation run the h(x) describes the interface fluctuations. The latter can be Fourier-transformed leading to the following expression for power spectrum of the quasi-one dimensional where h q is the one-dimensional Fourier transform of h(x) with q as the wave number, denotes the time average, k B is the Boltzman constant, and L x ,L z are width and thickness of the simulation box. The interfacial stiffness, ∼ γ m , is used as a fair estimation of the γ m , whereas different orientations of the crystal structure enable the determination of γ m anisotropy. Nevertheless, at this point, it is worth mentioning that the direct studies of model 57,59,60,76 and realistic 62,67,77,78 systems suggest that this effect is usually relatively weak, and therefore γ m can be obtained from ∼ γ m determined from a single crystallographic orientation. Then, ln h q 2 is a linear function of ln q with a slope equal to − 2 57 , and the intercept is directly related to ∼ γ m . Thus, γ m can be estimated by fitting the obtained dependence of h q 2 on q (expressed in logarithmic scales) to the linear function with the constant slope equal to − 2 and analyzing its intercept, see Fig. 2b, where discussed fits are presented for both studied systems. As the initial stage of this analysis we determined the region at which h q 2 is characterized by the linear dependence on ln q with slope equal to − 2.  Fig. 3. It can be clearly seen that N and U possess maxima located close to each other for both studied systems. Interestingly, the above maxima are located around the temperature at which the given system crystallizes during cooling. Arrows in Fig. 3 indicate the mentioned temperatures. Thus, CNT fairly predicts the thermodynamic conditions of crystallization. However, at this point, it is worth noting that despite the fact that both systems are very similar and that the procedures of the performed experiments are identical the system II crystallizes at temperatures much lower than T m ( T m − T = 45K for system I, whereas for system II this difference equals 20K ). This observation is even more intriguing when one takes into account that around 20K below the melting temperature, N for system II is about 3 times higher than for system I. Additionally, at discussed thermodynamic conditions U for system II is also much higher than for system I, and therefore from the CNT point of view, neither N nor U suspend the crystallization. Consequently, the system II should easily crystallize much faster than it is observed during cooling experiments. Moreover, we would like to put the reader's attention on another intriguing fact. Examining the cooling procedure, one might observe that despite the fact that at T cr the system

Discussion
Nevertheless, it must be noted that the crystallization process starts from the stochastic formation of the critical nuclei within the liquid. Therefore, to examine the crystallization tendency in detail, and then to confirm that the characteristics of crystallization process for examining systems are indeed similar we simulate the liquid structure at temperature 5K higher than T cr for 200ns or till the time at which we observed the crystallization event. The results are presented in Fig. 4a, where one can see that increase in the temperature implies an extension of the time needed for registration of the crystallization for both systems. This observation corresponds with a prediction of the CNT. Additionally, it is worth mentioning that system I is more sensitive for applied temperature changes because 3 from 5 simulation runs did not end in the crystallization, whereas system II always crystallized. Summarizing, we can state that at temperature equals T cr + 5K the time needed for the registration of crystallization is slightly shorter for system II. However, at discussed temperature conditions the nucleation rate is more than 25 times higher for system II than for system I. It implies that assuming that the time needed for crystallization is equal to about 100ns for system II (see Fig. 4a), we could anticipate that system I would persist in liquid state for 2500ns (the total time would be 25 time longer than 100 ns). However, within only 200ns , the system I crystallized twice. Thus, comparing both systems the prediction of CNT does not correspond to the observed results. Subsequently, we simulated both systems at temperature 5K lower than T cr . The results are shown in Fig. 4b. As one can observe both systems always crystallize within 5 ns, i.e., within the time which is 2 times shorter than in the case of cooling experiment. Similar to previous results the crystallization process proceeds slightly faster for system II. However, in this case, the differences between both systems are less prominent. Hence, at temperatures 5K lower than T cr the stability behaviors of studied systems can be considered as comparable. Moreover, we would like to note that CNT predicts that at temperature lower than T cr the crystallization process for system II should slow down due to the decrease in N and U , which is not observed in the performed experiments, see Fig. 5b.
At this point it has to be noted that N values presented in Fig. 3 are expressed in the unit of 1/nm 3 s , which implies that N , and hence the number of the created critical nucleuses depends on the system size. The two studied systems, which are comprised of the same amount of the molecules, are simulated at various temperatures. Therefore, they exhibit different volumes. Nevertheless, as it can be seen in Fig. 4, the differences in υ equal only about 5% and therefore its impact on N can be neglected.  1) and (2), the possible explanations of observed differences between prediction of CNT and computational experiment can be found. Both equations consider the diffusion of the molecules. It immediately implies that systems possessing a higher value of D exhibit higher values of N and U . It is crucial in the case of comparison between systems which are characterized by the significant differences in T m because the diffusion strongly depends on the temperature. The higher T m implies the faster diffusion of the system and consequently the greater N and U values are predicted by CNT. The latter seems to be crucial especially in the case of the system characterized by similar structure. In Fig. 5a one can see that D(T m ) for system II is higher for about one decade than for system I.
The latter immediately implies 10 times higher values of N and U for system II. Interestingly, one can see in Fig. 3 that N and U between both systems differ also about 10 times. Hence, the reported variation in CNT predictions for both systems is consistent with differences in D . Following this observation, in Fig. 5b we present N/D values for both systems, which indeed are very similar. This finding confirms that the main reason for divergences in CNT predictions for studied herein system is the noticeable difference in D values. At this point it is worth mentioning that for liquid close to the melting conditions the ratio between rotational and translational diffusion is constant and independent on the temperature 80 . Hence, neither translation nor rotation can be treated as a limiting factor for crystallization process of system I.
Summarizing, in this paper we calculate the N and U curves according to CNT for two RM systems, which differ exclusively in the value of the dipole moment. The use of proposed model molecules enables entire elimination of the molecular structure role in the crystallization process. Importantly, it makes also that, in contrast to standard simple model systems, obtained results for the system I cannot be uses to reproduce the results determined for the system II. It is also worth noting that, we calculate γ using CFM, instead of estimation of its value. Our results show that N and U curves differ strongly for two studied system. The system with higher value of the dipole moment is characterized by about 10 times higher N and U . Interestingly, despite the fact that the system II exhibits drastically higher values of N and U , it does not crystallize at expected thermodynamic conditions i.e., at conditions at which N and U for second system are sufficient to observe the crystallization process. Our results suggest that the main reason for observed discrepancies between results of performed computational experiments and CNT predictions is the diffusion constant.

Methods
We employ the previously proposed the quasi-real molecules of the rhombus shape, i.e., rhombus-like molecules (RMs), which remarkable advantage is that keeping the simplicity of classical model systems, they display the structural anisotropy typical for the real molecules and simultaneously enable the creation of the differently oriented dipole moments, µ 51,81-83 . On the basis of the results reported in Ref. 51 , we know that only one of the 5 different systems, which vary in the values and orientation of µ , crystallizes. Therefore, we use this system as a reference one. It consists of 4 identical atoms (of carbon atom mass) arranged, as we already mentioned, in a rhombus shape, which implies that RM possess short and long molecular axes (along diagonals of rhombus) simultaneously keeping identical bonds lengths. The latter is set to equal 0.14982 nm, which is close to 0.14 nm, i.e., a length of the bond linking two carbon atoms in the benzene ring. Additionally, the angles between bonds in RM are established to make one diagonal two times longer than the other. To ensure the best mimic of the real molecules by RM, the stiffness of bonds, angles, and dihedrals, as well as the non-bonded interaction between atoms of different RM molecules, are defined by OPLS all-atom force field parameters 84 provided for carbon atoms of the benzene ring. Then, the permanent µ is obtained by redefining charges of given atoms, i.e., those arranged along the longer axis are set to 0.0e ( e is an elementary charge), whereas those places along shorter one www.nature.com/scientificreports/ equal ±0.5e . In this way, we obtain the reference system I (i.e., system C2 from Ref. 51 ). The second examined herein system, i.e., system II, is identical to the previous one except of the difference in charges' values, which are set to ±0.75e for system II. Consequently, the two model systems differ only in the value of the dipole moment, µ II = 3 2 µ I . Schemes of two studied molecules are presented in the insets of Fig. 1. At this point is worth recalling that consistently to our previous mention the applied method for RM molecules creation makes results obtained for one of them cannot be somehow used to reproduce the results for another one. However, the identity of the molecular structure is preserved.
The heating process was performed by the use of the GROMACS software [85][86][87][88] at conditions of constant temperature and pressure, which were controlled by the Nose-Hoover thermostat 89-91 and Martyna-Tuckerman-Tobias-Klein barostat 92,93 ( p = 1000bar ) respectively. The increase in the temperature between subsequent steps is equal to 5K . The first half of simulations, which last for 5 ns (the time step dt = 0.001 ps), was spared for equilibration of the system, whilst the data was collected during the second half.